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A  method  is  presented  for  calculating  potential  flows  in  infinite  channels.  Given  a  collection  of 
N  sources  in  the  channel  and  a  zero  normal  flow  boundary  condition,  the  method  requires  an  amount 
of  work  proportional  to  N  to  evaluate  the  induced  velocity  field  at  each  source  position.  Previous 
schemes  have  been  based  either  on  conformal  mapping,  which  experiences  numerical  difficulties 
with  the  domain  boundary,  or  direct  evaluation  of  the  Green’s  function.  Both  require  0{N J)  work. 
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1  Introduction 

The  evaluation  of  potential  fields  in  infinite  channels  arises  as  a  numerical  problem  in  several  areas, 
most  notably  electrostatics  and  fluid  dynamics.  The  governing  equation  is  the  Poisson  equation, 

A*  =  -£  (1) 

subject  to  an  appropriate  boundary  condition.  In  this  paper,  we  will  restrict  our  attention  to 
two-dimensional  models  and  will  consistently  use  the  terminology  of  fluid  dynamics.  In  viscous 
incompressible  flow,  the  left-hand  side  is  the  stream  function,  the  right-hand  side  is  the  vorticity, 
and  the  condition  imposed  on  the  boundary  is  that  of  zero  normal  flow 


where  the  velocity  field  u  is  given  by 


u  •  n  =  0  , 


/a* 

U  *  (,  dy  ’  dx  )  * 


In  terms  of  the  stream  function,  this  is  equivalent  to  specifying  homogeneous  Dirichlet  boundary 
conditions 

*  =  0  .  (4) 

We  will  concentrate  on  particle  models,  where  the  vorticity  field  is  discretized,  not  by  a  mesh, 
but  by  N  point  vortices, 

N 

S  *  •$(*-*<,»-»)  .  (5) 

«sl 

Here,  6  is  the  Dirac  ^-function  and  is  the  strength  of  the  ith  point  vortex  located  at  (x,,  yt).  In 
vortex  methods,  what  we  would  like  to  compute  is  the  stream  function  and/or  velocity  field  at  each 
particle  position.  In  the  absence  of  boundary  effects,  the  desired  results  can  be  obtained  from  the 
free-space  Green’s  function  for  the  Poisson  equation  (-^lnr)  as 

* In((ar,  -  x,)2  +  (y;  -  y ,-)*)  for  »  =  1,...,JV  ,  (6) 

i*  * 

^ (yi-yj,x:-Xi)  for  i  =  h  jN 


(  2t  (x,  -  *j)2  +  (y,  -  y,)2 


Note  that,  using  direct  summation,  the  calculations  (6)  and  (7)  require  an  amount  of  work  propor¬ 
tional  to  N 2.  To  overcome  this  obstacle,  a  variety  of  fast  “JV-body”  methods  have  been  proposed  in 
the  last  few  years,  which  reduce  the  computational  complexity  to  0(N\ogN )  or  O(N).  These  in¬ 
clude  particle-in-cell  methods  [l,-* 3],  astrophysical  tree  codes  [2,3],  series  expansion  methods  [15,16], 
and  the  fast  multipole  method  [5,9,10]. 

Remark  1.1:  It  is  dear  from  (6)  and  (7)  that  the  stream  function  and  velocity  field  awe  un-  - 

bounded  in  the  neighborhood  of  a  point  vortex.  In  [7],  Chorin  introduced  the  idea  of  replacing 

the  point  vortices  by  “vortex  blobs”  whose  induced  field  is  held  constant  within  a  small  neighbor-  Codes 

hood  of  the  source.  More  recent  work  by  Hald  [11],  Beale  and  Majda  [4]  and  others  has  shown  /or 
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that  higher  order  accuracy  can  be  obtained  by  using  different  approximations  for  the  local  field. 
Outside  a  finite-size  core,  however,  the  velocity  field  due  to  a  vortex  blob  in  all  of  these  methods  is 
simply  that  of  a  point  vortex.  Since  we  are  interested  in  reducing  the  computational  cost  of  vortex 
methods,  which  is  generally  dominated  by  far  field  interactions,  we  will  ignore  the  precise  nature 
of  the  local  interactions  and  will  continue  to  use  the  point  vortex  model. 

For  a  straight  channel,  the  fluid  velocity  cannot  be  obtained  as  in  (6)  and  (7).  The  ma'u 
difficulty  is  that  the  zero  normal  flow  condition  can  only  be  satisfied  by  an  infinite  image  system 
(Section  2),  making  direct  summation  over  a  collection  of  point  sources  impossible.  The  most 
commonly  used  technique  for  overcoming  this  problem  in  constrained  flows  is  that  of  conformal 
mapping.  By  converting  the  calculation  to  one  in  the  upper  half  plane,  the  boundary  condition  can 
be  imposed  with  one  image  per  particle,  and  the  potential  flow  computed  as  in  (6)  and  (7)  with 
only  double  the  number  of  point  vortices  (Fig.  1).  An  attractive  feature  of  this  approach  is  that 
the  fast  jV-body  algorithms  for  free-space  calculations  may  directly  be  applied. 


C 

A 


Figure  1 

Conformal  mapping  of  the  channel  to  the  upper  half  plane.  The  left-hand  limit  points  A  and  C  are 
mapped  to  the  origin  and  the  four  solid  vertical  line  segments  in  the  channel  are  mapped  to  the  four 
semicircles  in  the  upper  half  plane.  Two  representative  particles  are  marked  by  the  small  circle  and 
square.  The  sera  normal  flow  boundary  condition  is  easy  to  apply  with  the  method  of  images  (each  source 
is  simply  reflected  across  the  z-axis  and  given  opposite  strength).  Unfortunately,  there  is  much  stretching 
and  contraction  of  the  physical  domain  which  can  cause  practical  difficulties. 

There  are  two  objections  to  this  mathematically  reasonable  procedure.  In  a  channel  with  zero 
normal  flow  boundary  conditions,  the  velocity  field  induced  by  a  point  source  decays  exponentially 
along  the  length  of  the  channel.  But  the  free-space  Green’s  function  used  in  the  upper  half  plane 
decays  as  £,  losing  much  of  the  physical  behavior  of  the  solution.  In  fact,  the  physical  behavior 
is  expressed  by  the  mapping  itself,  which  for  the  strip  {0  <  y  <  *•}  is  simply  ex .  The  second 
objection  is  that  a  discretization  of  the  boundary  is  often  required  (e.g.  for  vorticity  generation). 
Conformal  mapping,  however,  is  well  known  to  experience  numerical  difficulty  when  the  derivative 
of  the  map  has  a  great  dynamic  range  [14,17].  This  is  clearly  observed  in  Fig.  1,  where  the  images  of 
equispaced  points  along  the  top  and  bottom  of  the  channel  are  points  whose  separation  is  growing 
(or  contracting)  exponentially.  It  would,  on  both  counts,  be  much  preferable  to  remain  in  the 
channel  itself.  To  do  this  we  will  first  need  to  replace  the  infinite  image  system  by  an  analytic 
expression  for  the  Green’s  function.  This  can  be  obtained  through  elliptic  function  theory.  In 
[6],  Choi  and  Humphrey  derive  expressions  for  both  the  infinite  channel  and  a  closed  rectangular 


domain.  With  this  expression,  the  velocity  held  can  be  obtained  in  a  manner  analogous  to  the 
IV-body  calculation  of  equation  (7).  Direct  summation,  of  course,  will  require  0{N2)  work. 

In  this  paper,  we  propose  a  new  algorithm  for  two-dimensional  potential  flow  in  infinite  channels. 
It  is  based  on  the  analytically  derived  Green’s  function,  and  requires  an  amount  of  work  proportional 
to  N  to  evaluate  all  pairwise  interactions. 


2  Green's  function  for  an  infinite  channel 


We  begin  by  developing  an  explicit  expression  for  the  velocity  field  induced  by  a  point  vortex  in 
an  infinite  channel.  The  domain  is  defined  to  be  the  strip  {0  <  y  <  H}.  We  refer  to  the  direction 
x  increasing  as  downstream  and  to  the  direction  x  decreasing  as  upstream.  We  will  use  complex 
notation,  equating  the  points  (z«,yi)  with  the  complex  numbers  Z{.  If  we  define  u  by 


U(2f) 


=  _L_ 

rr?  2»  a*  -  Zi  ’ 


(8) 


then 

u(a\,yi)  =  (— Im(u(zi)),Re(u(z,)))  (9) 

is  the  velocity  field  induced  by  a  collection  of  point  sources  with  strength  (j  located  at  the  points 
Zj  =  (xj,yj).  In  the  remainder  of  this  paper,  we  consider  the  calculation  of  u  rather  than  u  and 
will  abuse  notation  by  referring  to  u  as  the  velocity  field. 

Let  us  now  suppose  that  a  source  of  unit  strength  is  located  inside  the  channel  at  zq  and  that  z 
is  a  second  point  inside  the  channel  with  z  jL  zq.  In  order  to  satisfy  the  zero  normal  flow  condition 
along  the  top  and  bottom  of  the  channel,  we  introduce  the  infinite  image  system  shown  in  Fig.  2. 

Let  us  first  add  up  the  contributions  from  the  images  with  positive  strength,  located  at  zo+2jHi 
(j  =  — oo,...,oo).  The  velocity  field  Uj(z)  induced  by  these  images  is  given  by  the  expression 


«i(*)  = 


But  from  ([8],  p.  36)  we  have 


oo  . 

y*  1 

jh U  *-*<>  +  2iHi 

1  y,  1  (  _ 

z  -  zo  +  ^  z  -  zo  +  2  jHi  z 

1  ,  V''  2 (z  —  Zq) 

z-  Zo  "  (z  -  zo)3  +  4 B2j2 


1 

zo  —  2 jHi 


so  that 
where 


w  x  1  2z  ^  i 
coth(«)=-+T|:-rnEf, 


iii(z)  =  <r  •  coth(o(z  -  zq))  , 


<T  = 


2  H 


(10) 

(ID 

(12) 

(13) 

(14) 

(15) 
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zq  +  4Hi 


Zo  +  2  H  i 


Zq  +  2  Hi 


zo  —  2Hi 
So  —  2Hx 


Figure  2 

Enforcing  boundary  conditions  by  the  method  of  images.  Successive  reflection  across  the  top  and  bottom 
boundaries  creates  the  image  system  shown.  The  images  at  positions  Sq  +  2jHi,  j  ss  —  oo, oo  have 
strengths  of  opposite  sign  from  those  at  positions  Zq  +  2jHi. 


For  the  images  with  negative  strength,  located  at  So  +  2 jHi  (j  =  -oo, ...,  oo),  we  obtain  an  induced 
velocity  field  iij,  given  by 

U2(z)  =  — <r  •  coth(<7(z  —  z0))  .  (16) 

The  net  velocity  field  is,  therefore, 


u(z)  =  a  •  (coth(<r(z  -  zo))  -  coth(<7(z  -  20)))  • 

A  simple  integration  yields  the  stream  function  $  induced  by  a  point  vortex, 


A  different  derivation  of ¥  is  given  by  Choi  and  Humphrey  in  [6].  As  mentioned  previously,  with 
this  analytic  expression  for  the  pairwise  interaction,  the  evaluation  of  the  velocity  field  at  each  of 
the  N  source  positions  can  be  carried  out  in  0(N 2)  operations.  In  order  to  develop  a  fast  algorithm 
tailored  to  the  channel  itself,  we  need  to  examine  the  properties  of  the  Green’s  function  in  more 
detail. 


2.1  Upstream  and  Downstream  Expansions 

Let  us  suppose  that  z  is  downstream  of  zq  (Re(z  -  zq)  >  0).  Then 


coth(<7(z  —  Zq ))  = 


c<r  (*-*o)  c-<T-(*-*o) 


c<t  (x-*o)  _  e-<r-{x-*o) 
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,fc  T.fj  .1  : 


=  "I  + 


1  _  e-2»  (*-*o) 


=  -1  +  2-£< 


,2  a  zc  k  —2 oxk 


=  1  +  2  •  £  e2"10*  *  e~2°tk  •  (22) 

fc= l 

Note  that  (21)  can  be  obtained  from  (20)  only  if  e-2<T'(*-*°)  <  1  which  is  ensured  by  the  condition 
that  z  be  downstream  of  the  source.  From  (17),  then,  the  velocity  field  downstream  of  a  unit  source 
at  zq  has  the  expansion  (about  the  origin) 

u(z)  =  a  •  (coth(<r(z  —  zo))  -  coth(<7(z  -  zo)))  (23) 

=  2<r  •  f^(e29t°k  -  e2aS°k)  •  e .  (24) 

From  this,  it  is  immediately  obvious  that  the  decay  in  the  field  is  exponential  along  the  length  of 
the  channel.  The  main  reason  for  developing  an  expansion  of  this  form,  however,  is  that  it  allows 
us  to  effectively  use  the  superposition  principle.  By  this  we  mean  the  following: 

Theorem  2.1  Suppose  that  m  sources  with  strengths  {q},  j  =  l,...,m}  are  located  at  points 
{zj,  j  =  l,...,m},  with  Re(zj)  <  r.  Then  for  any  point  z  further  downstream  (Re(z)  >  r), 
the  velocity  u(z)  induced  by  the  sources  is  given  by 


where 


«(z)  =  YL ak  •  e~ 


=  2 v'22qj-(e2<r'’k  -e2'*’1')  . 


The  error  in  truncating  the  expansion  (25)  after  p  terms  has  the  bound 


u(z)  -£flre' 


A  •  zr+1 
~  1  —  * 


where 


A  =  4<t  •  £  \qj\  and  x  =  . 

i-i 


Proof:  The  coefficients  a*  are  obtained  directly  from  equation  (24).  The  error  bound  is  a  conse¬ 
quence  of  the  triangle  inequality  and  the  fact  that  (24)  expands  the  field  due  to  a  single  source  as 
the  sum  of  two  geometric  series.  □ 

The  upstream  direction  is  treated  in  an  analogous  fashion.  If  Re(z  —  zo)  <  0,  then  the  velocity 
due  to  a  source  at  zq  can  be  expressed  as 

u(z)  =  2<t  •  £(e-2‘Tlofc  -  e-2<rtok)  •  e2'*k  .  (29) 


L 


Theorem  2.2  Suppose  that  m  sources  with  strengths  {qj,  j  =  1, m}  are  located  at  points 
{zj,  j  =  with  Re(rj)  >  r.  Then  for  any  point  z  further  upstream  (Re(z)  <  r),  the 

velocity  u (z)  induced  by  the  sources  is  given  by 


u  (z)-Zbt.e** 

fc= 1 

where 

bk  =  2o  ■  f;  qj  ■  -  e-'”*’k)  . 

j= i 

jTAe  error  in  truncating  the  expansion  (SO)  after  p  terms  has  the  bound 


**i 


A  •  iP+1 


1  —  x 


where 


A  =  4<t  •  [fy  l  and  x  as  . 


(30) 


(31) 


(32) 


(33) 


Definition  2.1  The  expansions  given  by  (SO)  and  (25)  will  be  referred  to  as  upstream  and  down¬ 
stream  expansions,  respectively.  For  a  given  collection  of  sources,  the  pair  will  be  referred  to  as 
5-expansions. 


The  representation  of  the  velocity  field  by  means  of  these  expansions  may  be  viewed  as  an 
analog  of  the  multipole  decomposition  of  the  field  due  to  a  collection  of  sources  in  free  space.  It  is 
important  to  keep  in  mind,  however,  that  both  their  rate  of  decay  and  region  of  convergence  are 
quite  different. 

Before  examining  the  properties  of  5-expansions  any  further,  we  demonstrate  their  usefulness 
in  computing  far  field  interactions  with  a  simple  example.  For  this,  suppose  that  U\  and  U?  are 
two  sets,  each  containing  N  point  vortices,  located  inside  a  channel  of  width  H,  and  separated  by  a 
distance  d  (Fig.  3).  To  compute  the  velocity  at  each  position  in  U\  due  to  the  sources  in  Z72  (or  the 


Figure  3 

Two  duster*  of  point  vorticea  located  inside  a  channel  of  width  H.  The  distance  between  the  two  dusters 
is  denoted  by  d. 


velocity  at  each  position  in  U2  due  to  the  sources  in  U\)  by  means  of  the  Green’s  function  would 
require  0(N2)  operations.  Let  us  instead  form  a  downstream  expansion  due  to  the  sources  in  U\ 


and  an  upstream  expansion  due  to  the  sources  in  U j.  From  (27)  and  (32),  it  is  easy  to  determine 
a  priori  how  many  terms  are  needed  to  achieve  a  relative  precision  of  e.  We  simply  require  that 
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p+ 1 


«  € 


or 


P  * 


H- Ini 
t  •  d 


(34) 


which  is  independent  of  N.  The  cost  of  formation  of  the  two  expansions  is  clearly  proportional  to 
Np.  Evaluating  the  two  expansions  at  all  points  in  the  relevant  cluster  also  requires  an  amount  of 
work  proportional  to  Np,  so  that  the  total  computation  scales  linearly  with  N,  assuming  that  the 
relative  precision  c  and  separation  distance  d  are  fixed. 


3  The  Shifting  Lemma  for  5-expansions 

The  fast  algorithm  to  be  described  depends  not  only  on  the  formation  and  evaluation  of  5- 
expansions,  but  on  their  analytic  manipulation.  The  following  obvious  lemma  describes  the  neces¬ 
sary  tools. 

Lemma  3.1  Suppose  that 

u(*u)=f (35) 

fcat  1 

and 

u(zd)  =  f^ak-e-2^k  (36) 

Arsl 

are  the  up  and  downstream  expansions  about  the  origin  due  to  m  sources  with  strengths  {qj,  j  = 
1, ...,  to}  which  are  located  at  points  {zj,  j  =  1,..., m},  with  — r  <  Re(2,)  <  r  for  some  r  >  0.  Then 

=  (37) 

fc= l 

and 

=  (38) 

*= l 

are  the  corresponding  up  and  downstream  expansions  about  zq,  where 

0k  =  bk  •  e3<rt°k  (39) 

and 

<*k  =  a/c  •  e~2<rx°k  .  (40) 

Furthermore,  the  error  bounds  for  the  shifted  S  -expansions  are  exactly  the  same  as  those  for  the 
original  S -expansions. 


Note  that  the  behavior  of  shifted  5-expansions  contrasts  sharply  with  that  of  multipole  expan¬ 
sions  in  free-space  (see  [9,10]).  In  the  latter  situation,  the  'validity  and  accuracy  of  an  expansion 
depends  not  only  on  the  source  positions  but  on  the  location  of  the  expansion  center.  Note  also 


from  (26)  and  (31)  that  the  coefficients  of  5-expansions  about  the  origin  are  pure  imaginary.  If  the 
centers  of  the  shifted  expansions  are  chosen  to  lie  along  the  z-axis,  then  the  coefficients  in  (39)  and 
(40)  are  also  pure  imaginary,  yielding  a  savings  in  both  computational  cost  and  storage. 

Remark  3.1:  To  this  point,  we  have  been  viewing  5-expansions  as  representations  of  the  far 
field  due  to  a  distribution  of  sources.  It  is  possible,  however,  to  view  them  in  a  different  light. 
The  expansions  (35)  and  (36)  of  the  preceding  lemma  are  valid  outside  the  strip  — r  <  Re(zj)  <  r. 
By  choosing  a  point  20  upstream  of  the  strip  boundary  (Re(zo)  <  —  r),  the  shifted  expansion  (37) 
yields  a  representation  of  the  induced  field  in  a  neighborhood  of  zq.  The  same  obviously  holds  for 
shifting  a  downstream  expansion  in  the  downstream  direction  (38).  These  are  local  representations 
of  the  field,  the  analogs  of  Taylor  series  in  free-space,  just  as  the  far  field  5-expansions  are  the 
analogs  of  multi  pole  expansions. 

4  The  Channel  Decomposition  Algorithm 

In  this  section,  we  describe  the  first  part  of  the  fast  algorithm.  The  basic  idea  is  to  subdivide  the 
channel  into  vertical  strips  and  to  use  5-expansions  to  compute  far  field  interactions. 

H 


(«)  (*) 

Figure  4 

Decomposition  of  the  channel  into  “elementary  stripe.*  The  original  distribution  of  particles  is  shown  in 
(a).  In  (b),  a  finite  domain  containing  all  particles  has  been  subdivided  into  rectangular  regions,  each  of 
which  has  an  aspect  ratio  of  one  third.  5-expansions  can  be  used  to  compote  the  interactions  between 
particles  contained  in  non-neighboring  strips. 


The  “elementary”  strips  into  which  the  channel  is  refined  have  an  aspect  ratio  of  one  third. 
The  reason  for  not  subdividing  too  much  further  is  clear  from  equation  (34).  As  d  approaches  0, 
the  number  of  terms  required  to  achieve  a  fixed  precision  grows  arbitrarily  large.  If  we  stop  using 
expansions  when  d  =  H/ 3,  however,  then  the  number  of  terms  required  is  given  by 


P~ 


3  •  In  7 


L  <  In  -  . 
c 


(41) 


We,  somewhat  arbitrarily,  choose  to  stop  subdividing  the  channel  at  this  point.  We  will,  of 
course,  need  to  compute  the  interactions  within  an  elementary  strip  and  between  nearest  neighbor 
strips.  This  part  of  the  calculation  will  be  described  in  section  5.  It  relies  on  some  additional 
analysis  and  the  Fast  Multipole  Method  for  free-space  problems. 

The  remainder  of  this  section  is  devoted  to  a  description  of  the  channel  decomposition  algorithm. 
The  main  strategy  used  is  that  of  clustering  particles  at  a  variety  of  spatial  length  scales  and 
computing  distant  interactions  by  means  of  5-expansions.  We  begin  by  determining  the  locations 
of  the  extreme  upstream  and  downstream  particles.  The  corresponding  section  of  the  channel  is 


considered  to  be  the  computational  domain,  and  a  sufficient  number  of  elementary  strips  are  created 
to  cover  the  region  (Fig.  4). 

We  proceed  by  introducing  a  binary  tree  structure  which  groups  particles  at  coarser  and  coarser 
levels  (Figure  5).  Level  0  corresponds  to  the  finest  discretization  of  space  (the  elementary  strips), 
while  level  /  +  1  is  obtained  from  level  /  by  the  merger  of  two  strips.  The  resulting  strip  at  the 
higher  level  is  referred  to  as  the  parent,  while  the  two  strips  being  merged  are  referred  to  as  the 
children.  Two  strips  at  the  same  refinement  level  are  said  to  be  nearest  neighbors  if  they  share 
a  boundary,  otherwise  they  are  said  to  be  well-separated.  By  construction,  then,  the  minimum 
distance  between  well-separated  strips  is  Hj 3,  and  in  order  to  achieve  a  precision  of  (  in  computing 
interactions  via  5-expansions  we  need  only  choose  the  number  of  terms  to  be  p  =  fln(l/e)].  At 
coarser  levels,  the  number  of  terms  can  obviously  be  decreased. 


(«)  (») 


(c)  (d) 

Figure  5 

In  (a.),  eight  elementary  strips  are  shown  which  cover  the  computational  domain.  This  level  of  spatial 
refinement  is  referred  to  as  level  0.  In  (b),  (c)  and  (d),  pairs  of  strips  are  successively  merged  to  form 
coarser  and  coarser  subdivisions  of  the  channel.  The  "center’’  of  a  strip  is  defined  to  be  the  midpoint  of 
the  segment  of  the  X-axis  bounded  by  that  strip,  as  indicated  in  (d). 

Definition  4.1  The  center  of  a  strip  is  defined  to  be  the  midpoint  of  the  segment  of  the  x-axis 
bounded  by  that  strip  (Fig.  5  (d)). 

Other  notation  used  in  the  description  of  the  algorithm  includes 

F/*t  a  p-term  upstream  expansion  about  the  center  of  strip  t  at  level  I,  describing  the  far  field  due 
to  the  particles  contained  inside  the  strip. 

Fft  a  p-term  downstream  expansion  about  the  center  of  strip  i  at  level  /,  describing  the  far  field 
due  to  the  particles  contained  inside  the  strip. 

a  p-term  local  5-expansion  (see  Remark  3.1)  about  the  center  of  strip  i  at  level  /,  describing 
the  field  due  to  all  particles  upstream  of  strip  i’s  nearest  neighbors. 

Lft  a  p-term  local  5-expansion  (see  Remark  3.1)  about  the  center  of  strip  i  at  level  /,  describing 
the  field  due  to  all  particles  downstream  of  strip  i’s  nearest  neighbors. 


Interaction  list  for  strip  i  at  level  /,  it  is  the  set  of  strips  which  are  children  of  the  nearest 
neighbors  of  i’s  parent  and  which  are  well-separated  from  strip  i  (Figure  6). 


Figure  6 

The  interaction  list  for  strip  i  at  level  1.  Strips  marked  with  a  “u”  are  upstream  members  of  the  list,  while 
those  marked  with  a  “d”  are  downstream  members  of  the  list.  Note  that  thick  lines  correspond  to  mesh 
level  1+1. 

The  channel  decomposition  algorithm  is  a  two-pass  procedure.  In  the  first  (upward)  pass,  we 
form  the  far  field  5-expansions  and  Ffi  for  all  strips  at  all  levels,  beginning  at  the  level  of 
elementary  strips.  In  the  second  (downward)  pass,  we  form  the  local  5-expansions  Lft  and  Lfi  for 
all  strips  at  all  levels,  beginning  at  the  coarsest  level. 

To  see  how  the  latter  part  is  accomplished,  suppose  that  at  level  l  +  1,  the  local  expansions  Lu 
and  Ld  have  been  obtained  for  each  strip  i.  Then,  by  using  lemma  3.1  to  shift  these  expansions 
to  the  centers  of  strip  i’s  children,  we  obtain  up  and  downstream  expansions  for  each  child  strip 
at  level  l,  describing  the  velocity  field  due  to  all  particles  up  and  downstream  of  strip  i’s  nearest 
neighbors.  For  each  strip  j  at  level  l,  then,  the  interaction  list  is  precisely  that  set  of  strips  whose 
contribution  to  the  potential  must  be  added  in  order  to  create  Lfj  and  Lfj  (Fig.  6).  For  each 
upstream  member  of  the  list,  we  use  Lemma  3.1  to  shift  the  center  of  the  corresponding  far  field 
expansion  Fu  to  the  center  of  strip  j  and  add  the  result  to  the  upstream  expansion  obtained  from 
the  parent.  Similarly,  for  each  downstream  member  of  the  list,  we  use  Lemma  3.1  to  shift  the 
center  of  the  corresponding  far  field  expansion  Fd  to  the  center  of  strip  j  and  add  the  result  to  the 
downstream  expansion  obtained  from  the  parent.  Note  that  at  the  coarsest  level,  Ld  and  Lu  are 
equal  to  zero,  since  there  are  no  well-separated  strips  to  consider. 

Finally,  for  each  strip  j  at  the  finest  level,  we  evaluate  the  local  expansions  Ld:  and  Lqj  at  the 
position  of  each  particle  contained  in  the  strip. 

Algorithm  1 

Comment  [Set  number  of  terms  to  be  used  in  expansions.] 

Choose  the  precision  t  to  be  achieved.  Set  the  number  of  terms 

in  all  expansions  to  p  =  fln(l/«)] . 

Upward  Pass 

Step  1. 

Comment  [Decompose  the  channel  into  elementary  strips.] 


Define  elementary  strip  width  to  be  Swu  =  if/3. 

Compute  xmin  =  x-coordinate  of  extreme  upstream  particle  position. 
Compute  xmax  =  x-coordinate  of  extreme  downstream  particle  position. 
Compute  number  of  elementary  strips  K  =  f(rma»  —  xm<n)/S«,,Y|. 
Compute  height  of  binary  tree  nlev  =  flog2  K\ . 

Step  2. 

Comment  [Form  far  field  5-expansions  at  finest  level.] 
do  i  =  1, ...,  K 

Form  p-term  up  and  downstream  expansions  F^i  and  Fd{ 
by  using  Theorems  2.1  and  2.2. 
end  do 

Step  3. 

Comment  [Form  far  field  5-expansions  at  all  coarser  refinement  levels.] 
do  /  =  1,  ...,n/e v 

Form  p-term  up  and  downstream  expansions  end  Fdi 
for  each  strip  i  at  level  l  by  using  Lemma  3.1 
to  shift  the  center  of  each  child  strip’s  expansions  to  the  current 
strip  center  and  adding  them  together, 
end  do 

Downward  Pass 


Step  4. 

Comment  [Form  local  5-expansions  at  all  refinement  levels.  Recall  that  Lu  and  Ld  are  zero  at  level  nlev 
since  there  are  no  well  separated  strips  to  consider.] 

do  l  =  nlev  —  1,  ...,0 

For  each  strip  i  at  level  I,  initialize  L“<  and  Ld{ 
by  shifting  the  Lu  and  Ld  expansions  of  strip  i’s  patent  to  the 
center  of  strip  t.  For  each  strip  in  i’s  interaction  list,  determine 
whether  it  is  up  or  downstream  of  strip  i.  If  upstream,  shift  the 
center  of  the  corresponding  Fu  expansion  to  i’s  center  and  add  to 
Ly{.  If  downstream,  shift  the  center  of  the  corresponding 
F ■  expansion  to  i’s  center  and  add  to  Lf(  (Fig.  6). 
end  do 

Step  5. 

Comment  [Local  5-expansions  are  now  available  at  the  finest  mesh  level.  They  can  be  used  to  compute 
the  velocity  field  due  to  all  particles  outside  the  nearest  neighbor  elementary  strips.] 

do  i  =  1, ...,  K 

For  each  particle  located  in  elementary  strip  i,  evaluate 
Lq  {  and  Lq  {.  Add  results  together, 
end  do 


A  brief  operation  count  of  the  channel  decomposition  algorithm  follows. 


Step  Number  Operation  Count  Explanation 


Step  1 

order  N 

examine  each  particle  position  to  determine  extreme 
up  and  downstream  coordinates. 

Step  2 

order  2  Np 

Each  particle  contributes  to  an  upstream  and 
a  downstream  expansion. 

Step  3 

order  K  p 

The  number  of  nodes  in  a  binary  tree  is  less 
than  twice  the  number  of  leaves,  so  that  the  total 
number  of  nodes  is  of  the  order  K .  For  each  node, 
an  amount  of  work  of  the  order  p  is  performed. 

Step  4 

order  5 p  ■  K 

For  each  strip  at  each  level,  there  are  at  most  three 
entries  in  the  interaction  list.  For  each  entry,  the 
amount  of  work  is  proportional  to  p.  In  addition, 
two  p-term  expansions  must  be  obtained  from  the  parent. 

Step  5 

order  2  Np 

Two  p-term  5-expansions  are  evaluated 
for  each  particle. 

The  estimate  for  the  running  time  is  therefore 

N  ■  (4p+  1)  +  K  6p.  (42) 

5  The  Evaluation  of  Nearest  Neighbor  Interactions 

The  channel  decomposition  algorithm  has  left  us  with  a  sequence  of  uncoupled  problems  to  consider. 
For  each  elementary  strip,  we  must  compute  the  internal  interactions  as  well  as  the  effects  of  the 
sources  contained  in  that  strip  on  the  particles  in  the  nearest  neighbors. 


Figure  7 

In  the  second  put  of  the  algorithm,  interactions  ue  computed  within  each  elementary  strip  and  between 
neuest  neighbors.  This  is  accomplished  by  muching  along  the  channel,  considering  one  strip  at  a  time, 
and  accounting  for  its  influence  on  all  relevant  puticles. 

Because  of  their  poor  convergence  rates  in  this  regime,  5-expansions  are  of  limited  use.  We 
could  proceed  by  direct  evaluation  of  the  remaining  interactions  through  the  use  of  the  Green’s 
function,  but  the  asymptotic  complexity  of  such  an  algorithm  would  be  0(N2).  Let  us  instead 
examine  one  of  the  subproblems  in  more  detail. 


We  begin  by  reconsidering  the  method  of  images  used  to  impose  the  zero  normal  flow  condition 
in  Figure  2.  Successive  reflection  across  the  top  and  bottom  of  the  channel  yields  a  one-dimensional 
array  of  squares  (Fig.  8).  These  are  either  copies  of  the  channel  section  itself  or  of  its  reflection 
across  the  bottom  boundary,  offset  by  2 jHt  for  some  integer  j.  Note  that  we  are  only  acting  on 
the  sources  contained  within  the  central  elementary  strip,  but  that  we  will  compute  the  velocity 
field  at  particle  positions  within  all  three  elementary  strips  of  which  the  square  is  composed.  In 
this  manner,  all  interactions  will  have  been  accounted  for  exactly  once. 

The  problem,  again,  is  how  to  account  for  the  sources  in  all  image  squares.  We  present  a 
solution  based  on  multipole  expansions. 


H 


Figure  8 

The  channel  section  and  its  translated  images  are  represented  by  boxes  labelled  C.  The  square  obtained 
by  reflection  across  the  bottom  boundary  and  its  translates  are  labelled  C. 


5.1  Multipole  Expansions 

We  will  require  two  results.  For  the  first,  suppose  that  m  point  vortices  with  strengths  g,  and 
positions  z,-  are  located  within  a  disk  of  radius  r  centered  at  the  origin.  Then,  for  a  point  z  with 
|z|  >  r,  the  velocity  field  v(z)  induced  by  the  sources  is  given  by  a  multipole  expansion  of  the  form 

v(*)  =  £  ff  .  (43) 

where 

«*  =  -**"1  • 

1=1 

The  error  in  truncating  the  sum  after  s  terms  is 


(44) 


where 


(46) 


^  =  23  Ittl  and  C  = 

i=l 

For  a  proof,  see  [9]. 

Note  that  in  order  to  obtain  a  relative  precision  of  e  (with  respect  to  the  total  charge),  the 
number  of  terms  required  in  the  series  representation  of  v  is  approximately  —  logc(e),  independent 
of  m,  the  number  of  source  charges. 

The  second  result  we  need  is  contained  in  the  following  lemma,  which  describes  the  conversion 
of  a  multipole  expansion  into  a  local  (Taylor)  expansion  inside  a  circular  region  of  analyticity. 


Lemma  5.1  (Conversion  of  a  Multipole  Expansion  into  a  Local  Expansion)  Suppose  that 
m  sources  of  strengths  qi,qi, ...,  qm  are  located  inside  the  circle  D\  with  radius  R  and  center  at  z0, 
and  that  |zo|  >  (c  +  1)-R  with  c  >  1.  Then  the  corresponding  multipole  expansion 

vW=£<^j1,  (47) 

converges  inside  the  circle  of  radius  R  centered  about  the  origin.  Inside  Z)2»  the  potential  due 
to  the  charges  is  described  by  a  power  series: 


j=o 


(48) 


where 


Furthermore,  for  any  s  >  max  ^2,  ,  an  error  bound  for  the  truncated  series  is  given  by 

v4(4e(s  +  c)(c  +  1)  +  c3) 


v(*)-£V*' 

1=0 


c(c  -  1) 


(50) 


where  A  is  defined  in  (46)  and  e  is  the  base  of  natural  logarithms. 

Proof:  See  [9],  □ 


Definition  5.1  Two  squares  with  sides  of  length  2d  are  said  to  be  well-separated  if  they  are  sepa¬ 
rated  by  a  distance  2d. 

Remark  5.1:  Let  A  and  B  be  well-separated  squares  with  sides  of  length  2d,  and  let  Da  and 
Db  be  the  smallest  disks  containing  the  boxes  A  and  B ,  respectively.  Then  the  disks  have  radii 
y/2  •  d,  and  the  distance  from  the  center  of  one  disk  to  the  closest  point  in  the  other  disk  is  at  least 
(4  —  y/2)  -d.  Letting  c  =  (4  -  \/2)/v/2  ss  1.828,  the  error  bound  (50)  applies  with  a  truncation  error 
using  s-term  expansions  of  the  order  c“*. 

Remark  5.2:  In  this  section,  the  center  of  a  square  refers  to  its  geometric  center  and  not  to  its 
strip  center  (Definition  4.1). 


5.2  Reduction  to  a  Free  Space  Problem 

We  will  use  Lemma  5.1  to  account  for  all  image  sources  outside  the  nearest  neighbor  squares.  The 
remaining  calculation  can  then  be  carried  using  the  free  space  Green’s  function  (see  page  1).  We 
begin  by  choosing  a  coordinate  system  with  the  origin  lying  at  the  center  of  Co.  For  each  square 
Cj,  the  multipole  expansion  induced  by  the  contained  sources  is  of  the  form 

V(Z)  =  S(^F'  (51) 

where 

Zj  =  2 jHi  (52) 

is  the  square’s  center.  Note  that  the  coefficients  a*  of  such  a  multipole  expansion  are  translation 
invariant;  i.e.  they  are  identical  for  all  integer  j.  Moreover,  for  j  jt  0,  Cj  is  well-separated  from  Co, 
and  the  field  induced  inside  the  channel  is  accurately  representable  by  an  s-term  local  expansion, 
where  s  =  f— logc(c)]  is  the  number  of  terms  needed  to  achieve  a  relative  precision  «  (see  Remark 
5.1).  This  local  representation  is  given  by  Lemma  5.1  as 

*A*)  =  52  bm  '  (53) 

msO 

with 

*)(-*)*•  <«> 

Let  S  be  the  set  of  non-zero  integers.  To  account  for  the  field  due  to  all  well-separated  images 
Cj,  we  compute  the  coefficients  of  a  local  representation  by  adding  together  the  shifted  expansions 
of  the  form  (54)  for  all  Zj  with  j  €  S  to  obtain 

*(*)=  £Ca'*m  (55) 

mcO 

where 

**“  ■  r  1)(_i)‘  -  i56> 

The  summation  over  5  for  each  inverse  power  of  Zj  can  be  precomputed  and  stored.  For  powers 
greater  than  one,  the  series  is  absolutely  convergent.  For  (m  +  k)  =  1,  however,  the  series  is 
not  absolutely  convergent,  and  the  computed  value  depends  on  the  order  of  addition.  Choosing 
a  reasonable  value  for  the  sum  of  the  series  requires  consideration  of  the  physical  model.  For 
this,  suppose  that  the  only  particle  in  the  simulation  is  a  source  of  unit  strength  located  at  the 
origin.  Then  the  image  system  corresponds  to  a  uniform  one-dimensional  lattice,  and  by  symmetry 
considerations,  the  induced  velocity  at  any  lattice  point  must  be  zero.  But  the  net  velocity  of  the 
particle  at  the  origin  corresponds  to  the  summation  over  5  of  1  jzj,  so  that  we  set 


Lemma  5.2  For  k  >  1, 


E?r 

S  2 


Hk 


if  k  is  odd 


TjTrC(fc)  if  k  is  even 


(58) 


where  ((z)  is  the  Riemann  zeta  function. 
Proof:  For  k  odd,  we  have 


For  k  even,  we  have 


= 


s 


jr[(2jHi)k^  f^{-2jHi)k 
0  . 


OO  ^  OO  1 

1  yS,  1 

2k~1Hkik  ££  jk 

(-1  )*  y  1 

2*“*#*  jb  J*  ’ 


(59) 

(60) 

(61) 

(62) 

(63) 


To  account  for  the  well-separated  images  of  Co,  we  will  require  the  corresponding  multipole 
expansion.  It  is  easy  to  verify,  however,  that  for  such  squares,  centered  at  a  point  Wj,  the  expansion 
is  of  the  form 

v(z)  =  V]  7 — — rr  • 

^  (z  -  Wj)k 


(64) 


where 

7 k  =  -a*  •  (65) 

Except  for  (?o  and  <?j,  all  of  these  images  are  separated  from  Co,  and  as  above,  the  fields  they 
induce  inside  the  channel  section  are  accurately  representable  by  a  local  expansion, 


p 

£ 

mm  0 


Z6m-z 


with 


The  well-separated  images  C3  clearly  have  centers 

...,  —5 Hi,  —3Hi,3Hi,5Hi,  ... 
Let  T  be  the  set  of  integers  of  the  form 

{±(2.7  +  1),  j  =  1,2,  ...,oo) 


(66) 


(67) 


(68) 


(69) 


We  again  account  for  the  field  due  to  all  well-separated  images  by  forming  the  coefficients  of  a  local 
representation 


m=0 


where 


c total  _m 
°m  •  z 

y 

(70) 

(?<+*)■ 

(71) 

The  summation  over  T  for  (m  +  k)  >  1  is  absolutely  convergent.  For  (m  +  k)  =  1,  the  series  is 
not  absolutely  convergent,  but  symmetry  considerations  again  dictate  the  choice 


E— =  °- 

in 


T  ^ 


Lemma  5.3  For  k  >  1, 


1  =! 0 

Wj  |  12lz 


if  k  is  odd 


where  Bk  is  the  kth  Bernoulli  number. 
Proof:  For  k  odd,  we  have 


For  k  even,  we  have 


The  result  now  follows  from  the  equality  ([8],  p.  7) 


Y' _ l _ (22t  ~  1)y2*in  | 

~  (2j  +  l)2fc  ~  2  •  (21:)!  1  2fc| 


(72) 


(73) 


oo  ^  °°  i 

§  ((2j  +  1  )Hiy  +  2  (-(2 j  +  1)^0* 

(74) 

=  0  . 

(75) 

S* 

T  ) 

OO  .  OO  1 

-  Y  1  ,  V  1 

(76) 

%  ((2 j  +  \)HiY  fr[  (-(2 j  +  l)l?i)fc 

2  ^2,  i 

Hkik  i=[  (2i  +  *)* 

(77) 

2-(-l)*/J^  1 

Rk  ,tt(2i  + 1)** 

(78) 

(79) 
>*o 

□ 

If  we  add  the  computed  coefficients  flotai  from  (71)  to  the  coefficients  b)ftal  from  (56),  we  obtain 
a  single  local  expansion  which  describes  the  field  due  to  all  sources  outside  the  nearest  neighbor 
squares  of  Cq.  This  local  expansion  can  then  be  evaluated  at  all  particle  positions  in  Co- 
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The  Anal  step  in  the  algorithm  is  to  compute  the  velocity  field  due  to  the  free-space  sources 
within  Co,  Cq  and  Cj.  This  problem  is  handled  by  the  Fast  Multipole  Method  (FMM),  which 
requires  an  amount  of  work  proportional  to  n  +  m  to  evaluate  the  field  induced  by  n  sources  at  m 
points. 


Algorithm  2 


Comment  [Set  number  of  terms  to  be  used  in  expansions.] 

Choose  the  precision  e  to  be  achieved.  Set  the  number  of  terms 
in  all  expansions  to  s  =  floge(l/c)"|. 

Comment  [From  Algorithm  1,  we  are  given  that  the  number  of  elementary  strins  is  A'.] 

Define  n<  to  be  the  number  of  particles  in  the  if  A  strip. 

Clearly,  ni  +  n?  +  ...  +  n*  =  N,  the  total  number  of  particles. 


Comment  [Process  each  elementary  strip  .] 
do  «  =  1, ...,  K 

Define  C0  to  be  the  square  whose  central  third  is  strip  i. 

Step  1 

Form  coefficients  at  of  s-term  multipole  expansion  about  center  | 

of  C0  induced  by  sources  in  strip  i.  i 

Form  coefficients  yt  of  s-term  multipole  expansion  for  square  Co  ] 

via  equation  (65).  1 

Step  2  j 

Form  coefficients  6*  +  St  of  s-term  local  expansion  about  the 

center  of  Co  which  describes  the  field  induced  by  all  reflected  i 

sources  outside  the  nearest  neighbor  squares.  | 

Step  3 

Evaluate  local  expansion  at  all  particle  positions  in  strips  i  —  l,i  and  i  +  1.  I 

Step  4  j 

Compute  velocity  field  induced  by  sources  in  C0,  Co  and  Cj 
at  all  particle  positions  in  strips  i  —  1,  t  and  i  +  1  via  the  FMM. 
end  do  | 

A  brief  operation  count  of  the  Algorithm  2  follows. 

Step  Number  Operation  Count  Explanation 


Step  1  order  Nt  each  particle  contributes  to  an  s-term 

multipole  expansion  when  its  elementary  strip  is  being  processed. 

Step  2  order  Kt 2  The  creation  of  a  local  expansion  requires 

order  a2  work  and  is  carried  out  once  for  each  elementary  strip. 

Step  3  order  3 Nt  Three  local  expansions  are  evaluated  for  each  particle. 
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Step  4 


order  3a  N 


K  free-space  problems  are  solved,  each  of 
dimension  n«,  with  a  factor  of  three  included  to  account  for 
the  extra  image  sources  and  evaluation  locations.  The  factor 
a  represents  the  constant  for  the  linear  time  FMM. 

The  estimate  for  the  running  time  is  therefore 

N  •  (4s  +  3a)  +  K-s2.  (80) 

To  summarize,  then,  the  full  algorithm  consists  of 

1.  Decomposition  of  the  channel  into  elementary  strips, 

2.  Algorithm  1  to  compute  distant  interactions,  leaving  a  sequence  of  uncoupled  nearest  neigh¬ 
bor  problems  to  consider, 

3.  Algorithm  2  to  compute  nearest  neighbor  interactions. 


6  Numerical  Results 

A  computer  program  has  been  implemented  using  the  channel  decomposition  and  nearest  neighbor 
algorithms  of  this  paper.  For  testing  purposes,  we  randomly  assigned  particles  to  positions  within 
a  channel  section  of  length  5 H,  where  H  was  the  channel  width  (Fig.  4),  with  source  strengths 
between  0  and  1.  Five  digit  accuracy  was  requested  from  the  expansions.  In  the  first  part  of  the 
algorithm,  5-expansions  were  computed  to  10  terms,  while  in  the  second  part  of  the  algorithm, 
multipole  and  Taylor  expansions  were  computed  to  about  20  terms.  We  performed  the  calculations 
in  four  ways:  (1)  through  the  algorithm  of  this  paper  in  single  precision;  (2)  directly  from  the 
Green’s  function  in  single  precision;  (3)  directly  from  the  Green’s  function  in  double  precision; 
(4)  via  conformal  mapping  in  single  precision.  The  direct  evaluation  from  the  Green’s  function  in 
double  precision  was  used  as  a  standard  for  comparing  the  relative  accuracies  of  the  other  three 
methods  in  a  least  squares  sense.  Calculations  were  carried  out  on  a  SUN  3/50  workstation  using 
the  68881  co-processor. 

The  following  observations  can  be  made  from  Table  1. 

1.  The  accuracies  of  the  results  obtained  by  the  fast  algorithm  are  in  agreement  with  the  error 
bounds  given  in  this  paper.  In  fact,  the  results  are  consistently  more  accurate  than  either  of 
the  direct  calculations. 

2.  The  CPU  time  requirements  of  the  fast  algorithm  appear  to  grow  somewhat  superlinearly. 
The  reason  for  this  is  that  there  are  two  constants  asssociated  with  the  algorithm,  a  small 
one  for  the  channel  decomposition  and  a  larger  one  for  the  FMM.  The  observed  timings  are 
dominated  by  the  first  constant  for  100  and  400  particles,  and  by  the  second  constant  for 
the  larger  tests.  When  there  are  a  small  number  of  particles  per  strip,  the  FMM  with  its 
associated  overhead  is  simply  not  invoked. 


N 

*alg 

Tjir 

Tcm 

Ealg 

Edir 

Ecm 

100 

8.38 

34.8 

14.0 

4.5  10‘7 

7.2  10"7 

1.1  10- 

400 

53.1 

551 

223 

2.7  10-7 

4.1  10“7 

1.2  10- 

1600 

398 

(8820) 

(3550) 

4.3  10-7 

1.3  lO"6 

1.1  io- 

6400 

1890 

(141000) 

(56800) 

6.9  1Q“7 

5.2  lO*6 

3.4  lO" 

Table  1 

Table  of  CPU  times  in  seconds  required  by  the  fast  algorithm  (alg),  the  direct  Green’s  function  method 
(dir),  and  confoimil  mapping  with  direct  evaluation  of  the  resulting  jV-body  problem  (cm).  The  least 
squares  errors  for  the  three  methods  are  shown  in  the  last  three  columns.  Timings  in  parentheses  are 
estimated  by  computing  the  results  for  only  a  subset  of  100  of  the  particles.  The  corresponding  errors  are 
computed  from  that  smaller  data  set. 

3.  By  the  time  the  number  of  particles  reaches  6400,  the  fast  algorithm  is  about  75  times  more 
efficient  than  the  direct  Green’s  function  method. 

4.  Even  for  as  few  as  100  particles,  the  fast  algorithm  is  about  four  times  faster  than  the  direct 
calculation. 

7  Conclusions 

jA  fast  algorithm  for  potential  flow  in  channels  has  been  developed.  It  is  based  on  asymptotic 
expansions  which  we  refer  to  as  5-expansions,  some  analytic  observations  concerning  classical  mul¬ 
tipole  expansions  and  Taylor  series,  and  the  Fast  Multipole  Method.  The  asymptotic  CPU  time 
requirements  for  the  algorithm  grow  linearly  with  the  number  of  sources  and,  despite  its  complex 
structure,  numerical  experiments  demonstrate  that  dramatic  speedups  can  be  obtained  for  even 
moderate  size  particle  systems. 

In  its  current  form,  the  algorithm  requires  that  the  channel  boundaries  be  straight.  A  method 
applicable  to  channels  with  perturbed  boundaries  will  be  described  in  a  subsequent  paper.  '  * 
The  author  would  like  to  thank  V.  Rokhlin  for  several  useful  conversations. 


& 


I 


References 


'A  '  ■  P  hrr- 


[1]  C.  R.  Anderson,  A  Method  of  Local  Corrections  for  Computing  the  Velocity  Field  Due  to  a 
Distribution  of  Vortex  Blobs ,  J.  Comput.  Phys.,  62  (1986),  pp.  111-123. 

[2]  A.  W.  Appel,  An  Efficient  Program  for  Many-body  Simulation ,  Siam.  J.  Sci.  Stat.  Comput.,  6 
(1985),  pp.  85-103. 

[3]  J.  Barnes  and  P.  Hut,  A  Hierarchical  0(N  log  N)  Force-Calculation  Algorithm,  Nature,  324 
(1986),  pp.  446-449. 

[4]  J.  T.  Beale  and  A.  Majda,  Vortex  Methods  II:  Higher  Order  Accuracy  in  Two  and  Three 
Dimensions,  Math.  Comp.,  39  (1982),  pp.  28-52. 


Fa*,?  ,  /'•'  ./■;! 

if*. 


p 

$ 

ft* 

•  it* 


1 

% 

I 

■■iV 

.1 

li 


$ 

% 


1 


mJHUMiMiuin^inMJU,uuiiw?ww!wriiimnx7U'Aj'JUMmMTUi!uny  tuuinjiRUUAViuMMRUxuTfrmmv 


[5]  J.  Carrier,  L.  Greengard,  and  V.  Rokhlin,  A  Fast  Adaptive  Multipole  Algorithm  for  Particle 
Simulations,  Siam  J.  Sci.  Stat.  Comput.,  to  appear  (July,  1988). 

[6]  Y.  Choi  and  J.  A.  C.  Humphrey,  Analytical  Prediction  of  Two-Dimensional  Potential  Flow 
Due  to  Fixed  Vortices  in  a  Rectangular  Domain,  J.  Comput.  Phys.,  56  (1984),  pp.  15-27. 

[7]  A.  J.  Chorin,  Numerical  Study  of  S'ightly  Viscous  Flow,  J.  Fluid.  Mech.,  57  (1973),  pp.  785-796. 

[8]  I.  S.  Gradshteyn  and  I.  M.  Ryzhik,  Tables  of  Integrals,  Series,  and  Products,  Academic  Press, 
New  York,  1980  . 

[9]  L.  Greengard  and  V.  Rokhlin,  A  Fast  Algorithm  for  Particle  Simulations,  J.  Comput.  Phys., 
73  (1987),  pp.  325-348. 

[10]  L.  Greengard,  The  Rapid  Evaluation  of  Potential  Fields  in  Particle  Systems ,  MIT  Press,  Cam¬ 
bridge,  1988. 

[11]  O.  Hald,  Convergence  of  Vortex  Methods  for  Euler’s  Equations,  Siam  J.  Sci.  Stat.  Comput., 
16  (1979),  pp.  726-755. 

[12]  R.  W.  Hockney  and  J.  W.  Eastwood,  Computer  Simulation  Using  Particles,  McGraw-Hill, 
New  York,  1981  . 

[13]  A.  Leonard,  Vortex  Methods  for  Flow  Simulation ,  J.  Comput.  Phys.,  37  (1980),  pp.  289-335. 

[14]  S.  T.  O’Donnell  and  V.  Rokhlin,  A  Fast  Algorithm  for  the  Numerical  Evaluation  of  Conformal 
Mappings,  Technical  Report  554,  Yale  Computer  Science  Department,  1986. 

[15]  A.  M.  Odlyzko  and  A.  Schonhage,  Fast  Algorithms  for  Multiple  Evaluations  of  the  Riemann 
Zeta  Function ,  Trans.  Am.  Math.  Soc.,  to  appear  (July,  1988). 

[16]  L.  van  Dommelen  and  E.  A.  Rundensteiner,  Fast,  Adaptive  Summation  of  Point  Forces  in  the 
Two-Dimensional  Poisson  Equation,  submitted,  J.  Comput.  Phys. 

[17]  C.  Zemach,  A  Conformal  Map  Formula  for  Difficult  Cases,  J.  Comp.  Appl.  Math.,  14  (1984), 
pp.  207-215. 


21 


